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ABSTRACT 

We present a generic algorithm for generating Gaussian random initial 
conditions for cosmological simulations on periodic rectangular lattices. We 
show that imposing periodic boundary conditions on the real-space correlator 
and choosing initial conditions by convolving a white noise random field results 
in a significantly smaller error than the traditional procedure of using the power 
spectrum. This convolution picture produces exact correlation functions out to 
separations of L/2, where L is the box size, which is the maximum theoretically 
allowed. This method also produces tophat sphere fluctuations which are exact 
at radii R < L/4. It is equivalent to windowing the power spectrum with the 
simulation volume before discretizing, thus bypassing sparse sampling problems. 
The mean density perturbation in the volume is no longer constrained to be 
zero, allowing one to assemble a large simulation using a series of smaller ones. 
This is especially important for simulations of Lyman-a systems where small 
boxes with steep power spectra are routinely used. 

We also present an extension of this procedure which generates exact initial 
conditions for hierarchical grids at negligible cost. 



Subject headings: methods: numerical — large-scale structure of universe 



1. Introduction 

Cosmological N-body and hydrodynamic simulations have numerical resolution 
limitations imposed on both large scales comparable to the simulation box (Gelb and 
Bertschinger 1994, hereafter GB), as well as on small scales due to finite resolution (Hockney 
and Eastwood 1981). For simulations of Cold Dark Matter initial conditions, they found 
that for a box of width L=51.2 Mpc, the truncation of the power spectrum at the box 
length resulted in fluctuations on 8 Mpc spheres as which are systematically 40% lower 
than the analytical value. This factor of nearly 2 in error arises due to the truncation of 
the power spectrum while generating the initial conditions. They found that a box size 
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with L=100 Mpc was required to reduce the error to 10%. In cosmological applications, 
the actual value of ag is of significant interest since it corresponds to the mass fluctuations 
from which a rich cluster forms. In the Press-Schechter framework (Press and Schechter 
1974), cluster abundances will be statistically correct as long as <jr is exactly modeled 
in a simulation, where M = (4/3)nR 3 p c is the mass of the cluster in question. Since 
computational cost grows at least as fast as the third power of the simulation volume at 
fixed physical resolution, it is desirable to use numerical discretizations which allow the 
most accurate representation of physical processes with the minimal simulation volume. 
The first purpose of this letter is to demonstrate how the initial conditions in a periodic 
box can be chosen to eliminate this error entirely for tophat smoothing scales which are of 
size up to L/4. When an ensemble of simulations is performed, their average will reflect the 
correct tophat fluctuation distribution. 

The current trend in high-resolution simulations is to adaptively place high resolution 
grids in regions which are known to form high density objects (Bryan and Norman 1995, 
Dutta 1995, Navarro et al. 1995, Bartelmann and Steinmetz 1996). If this procedure is 
repeated recursively, one needs to be able to specify the initial conditions at very high 
resolution in these regions. The brute force approach needs to generate the initial conditions 
at high resolution everywhere, and throw out the regions which are never resolved. The 
convolution picture (Salmon 1996) allows a simple decomposition of the problem into 
multiple scales, each of which is rapidly computed using a FFT. This paper describes the 
optimal implementation of the convolution picture for cubical regular or hierarchical grids. 

A third related problem is the question of how one might embed a small periodic 
simulation volume inside of a larger one. For example, one might simulate a large lGpc 3 
volume of the universe at modest resolution, and identify a region in this simulation 
(perhaps an extremely rich cluster or empty void) which one would like to study at higher 
resolution. Since perturbations on scales of even 100 Mpc are small, one might envision 
drawing 100 Mpc box from the large volume and only simulating the smaller volume at 
higher resolution. If the simulation imposes periodic boundary conditions, the procedure 
described in this letter will allow us to generate such a box which agrees very well with 
the larger box in the central regions of overlap. The converse was studied by Tolman and 
Bertschinger (1996), which is how to run a large simulation at high resolution given a small 
one. This problem also fits in the framework of the convolution picture. 

The general solution to all three problems is to envision the generation of initial 
conditions in coordinate space instead of Fourier Space. We will initially lay down a 
random number at each grid point. Now gridpoints are uncorrelated with each other, 
and the correlator is a Dirac 5-function. Let us now convolve this initial condition with 
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some appropriate window W p (r). Convolution, being a linear process, will preserve the 
Gaussianity of the random field, but introduces a correlation function which is just the 
auto-correlator of W p . Numerically, this convolution will be performed again with the 
aid of FFTs. Since it is so similar to the Fourier method in implementation one can ask 
why it would give significantly different answers. The standard discretization of the power 
spectrum only samples at discrete frequencies k. For small k, this error can be large, 
especially if the slope of the power spectrum is steep. Discretizing a periodic correlator is 
equivalent to convolving the power spectrum with the window function of the simulation 
geometry, which smoothly averages over the power spectrum in the regions were it is poorly 
represented. For power spectra which have significant features, for example due to baryon 
oscillations, this is especially important. The correlator constructed in this manner should 
show absolutely no error for r < L/2. This is the qualitative reason for the improved 
performance of the convolution picture. 

The effect of truncating the power spectrum are particularly severe if the power 
spectrum is steep. As the slope n of the power spectrum P(k) = k n approaches n — > —3, 
most of the small scale top hat variance arises from large wavelength modes which are not 
included in the box size. For Cold Dark Matter power spectra, the slope approaches this 
value on small scales, and simulations of the early non-linear mass scales, for example in 
Lyman-a simulations, are more strongly affected by the loss of long waves (see for example 
Rauch et al. 1996 and references therein). 

This letter will mathematically describe the implementation in section 0. Numerical 
measures of the convolution performance are given in section |3] together with a discussion of 
its application to generate simulations which give the correct ensemble average. In section 
^ we describe how the convolution method allows the seamless adaptive generation of 
subpower for hierarchical grids, and how to generate consistent periodic sub- (super) boxes. 
Initial conditions on unstructured grids can be treated by the tree algorithm (Salmon 1996). 

2. Formalism 

We consider density fields 5(x) where x = (x,y,z) is a position in three dimensional 
Euclidean space. The autocorrelation function (or auto-correlator for short) is defined 
as £(r) =< 5(x)5(x + r) >. The density field can be transformed into Fourier space 
5(k) = V^ 1 J exp(— i\t-x)5(x.)d 3 x, where we measure the power spectrum P(k) = S(k)8(— k). 
The power spectrum is also the Fourier Transform of the auto-correlator. We define 
the spherical tophat window function Wt{u) = 3[sin(w) — ucos(u)]/u with which we 
measure the standard deviation of the density field smoothed on tophat spheres of radius 
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R a\ = 47r / P(k)WT{kR)k 2 dk. The traditional discretization to obtain a random density 
field is to choose a Gaussian random number for each mode 5(k) in Fourier space with 
variance P(|k|). For small k, the spacing of modes Ak ~ 2tt/L, where L is the simulation 
box length. Where k ~ 2tc/L, this sampling introduces significant coarse graining error, 
which we will call the truncation error. Within the Fourier domain, we next introduce 
the cubical tophat window W^k) = sm(2TrLk x ) sm(2TrLk y ) sin(2TTLk z )/[(2TrL) 3 k x kyk z ], 
which reflects the simulation geometry. To implement the convolution picture, we define 
the correlation kernel W p (k) = \J P{k). We call it a "picture" since it is equivalent to 
the Fourier space discretization for infinite volumes, so no new physics is introduced. 
Since P(k) is positive, so is W p , but not necessarily its inverse Fourier transform W p (x). 
One can verify that the auto-correlator is the auto-convolution of the correlation kernel 
£(r) =< W p (x)W p (x + r) >. We now generate a random field with zero correlation length 
<S (x) and define 5(x) = J 6 (x!)W p (x - x')d 3 x'. 

We now consider two alternatives to implement these relations on a periodic lattice. In 
Fourier space, periodicity means that only a discrete spectrum of modes are represented. 
At the modes which are represented, the traditional approach is to simply use the power 
spectrum P(k) at those discrete intervals. Two problems arise. If the power spectrum 
changes significantly between successive modes, one misses the intermediate power. The 
extreme case would be if the power spectrum had a feature, say a peak, which did not lie 
on one of the discrete grid points. The other problem is that the k = mode is now forced 
to be zero. We see that we really want some way of not just using the value of P{k), but 
some average over the width of the simulation volume We- 
ill the convolution picture, we impose periodicity on the correlator. We first set the 
correlator equal to zero outside of the box volume, and then create periodic images. We 
will add a superscript grid to indicate any quantity which has been made periodic using this 
convolution picture prescription. For the effective grid power spectrum, this is equivalent 
to replacing the simulation power by the full power spectrum convolved with Wc' 

pg" d (k) = J P(\k'\)W c (V - k)d 3 k'. (1) 

The effect of this convolution is biggest at small k, especially k = 0, where the analytic 
power spectrum is typically 0, but the windowed one takes on some finite value which is 
similar to the variance in the box simulation volume. The window function of the simulation 
geometry smoothly averages between the finite discrete sampling imposed by the periodic 
boundary conditions. Due to the anisotropy of Wc the grid power spectrum P gnd (k) has 
an angular dependence. We will discuss the effects of ([!]) in more detail in section || 
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3. Errors 

Figure [I] shows the ensemble averaged value of cr 8 computed using the traditional power 
spectrum as well as in the convolution picture. The GB estimate of cx 8 truncates the lower 
bound of the integral over the power spectrum in Equation (§) at k — 2n/L. The result for 
a model with power spectrum r = .25, cx 8 = 1 (Bardeen et al. 1998) is shown as the dotted 
line in the figure. In order to obtain the exact value of the periodic box cr 8 , we computed 
the actual tophat convolution on a discrete 128 3 grid without assuming a spherical cutoff in 
the power spectrum. This leads to the actual values of <j 8 shown as the dashed line. The 
solid line from the convolution picture was generated on a 256 3 lattice taken from a code 
which actually computes N-body simulation initial conditions. 

We notice the significant improvement in the linear theory analysis which is exact for 
box sizes greater than 32 h~ x Mpc. Even for box sizes of 16 hr x Mpc, where the 8 h~ x Mpc 
spheres start wrapping around the period volume, the convolution picture has an error of 
only 7%. By contrast, discretizing the power spectrum introduces an error of over 60%. We 
notice further that the convolution picture typically overproduces c g , while truncation in 
Fourier space underestimates it. We can understand it as follows. The correlation function 
is now forced to be periodic, so instead of steadily decaying at large separations, it increases 
beyond the periodicity radius. We have 

al = 4vr / £(r)W%(r)r 2 dr (2) 
Jo 

where W^(r) is the convolution of the tophat kernel with itself. It is zero for r > 2R, so we 
say that it has compact support. Periodicity thus causes (Q) to be overestimated because 
£ itself is overestimated for correlators which are positive and decreasing at r = L/2. 
Truncation in Fourier space, on the other hand, usually results in an underestimate because 
long waves are lost. 

A related problem is the box-to-box density variance. In the traditional power spectrum 
generation, the total density within each box is typically defined to be the cosmic mean 
density. In reality, this value varies even when averaged over a corresponding volume of the 
universe. The mean density of the box can be modeled as a change of the mean matter 
density O and incorporated in the evolution of the scale factor. This addresses one of the 
causes of systematically underestimating er 8 , which can result in systematic errors when 
simulating the abundance of virialized objects such as clusters of galaxies in finite volumes. 
In the convolution picture, the periodic correlator does not have zero mean, so the k = 
mode will have finite power. The variance thus induced will not be exactly equal to the true 
variance in cubical tophats of the same volume. The difference in shown in Figure 0. We 
see that the variance of the whole box variance is slightly overestimated, but typically only 
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Fig. 1. — <7g as a function of box size using various approximations. The dotted line was 
computed by GB94 using a truncated power spectrum. The dashed line is the actual tophat 
fluctuation in a 128 3 box. The solid line was obtained using the convolution picture. The 
latter is exact for box sizes greater than 32 h^ 1 Mpc. 
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by 10%. The tophat spheres of up to R = L/4 all have the exact variances, but the box 
variance itself corresponds to R ^ L/2. The correlator is already completely determined by 
requiring the smaller ones to be exact, so we have no additional freedom to adjust the large 
scale fluctuations. Luckily the corrections are small because the dominant contribution in 
(@) comes from R < L/4. 

A dual view of this picture is that the box-to-box variance is described by Equation 
([[]) where Wc is replaced by Wq. One could pose the question why one doesn't simply 
make such a replacement for all modes. The reason is that imposing exact variances in 
tophat spheres which can be embedded in half the simulation width, makes Equation ([!]) 
the unique solution. If one is interested in the dynamics and abundance of bound objects, 
the success of the Press-Schechter ansatz argues that the top-hat spheres should be of 
primary importance. A potential drawback of Equation ([TJ) is that the grid power may 
become negative, making it impossible to use it as a variance. On the numerical grid, one 
can add a constant to the numerical power spectrum to make it positive everywhere. This 
is equivalent to adding white noise, which shows up only in the r = bin of the two-point 
correlator. In practice, for the transfer function used in this study, the effect of varying 
£(0) between and £(Ar), the value of the correlator at one grid cell lag, did not affect 
measures like <jg by more than 1 part in 1000 even down to box sizes of 10 Mpc on a 
128 3 grid. We finally note that the simulation results must be analyzed in the same way 
as an observational survey. The power spectrum of the simulation volume is not a true 
representation of the infinite volume power spectrum. It is instead a windowed sample, 
which must be taken into account during statistical analysis. 



One direct advantage of the convolution picture is its adaptability to arbitrary domains. 
We first consider the situation where we want to add power on subgrid scales. Consider 
a fine rectangular grid (heavy lines) which is superposed on a coarser outer grid as shown 
in figure Also shown is the dotted extended fine grid, which is larger than the fine grid 
by an amount parameterized by R. We will assume that a random initial condition was 
already generated on the coarse grid. 

We first note that the correlation kernel can always be decomposed into a sum of two 
kernels W p = Wr + W' R where we can choose 



4. Subgrid Power 




if r > R 
if r < R 



(3) 
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r=0.25, cr R = 1 
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Fig. 2. — The standard deviation of the box density using the correlator (solid) and the 
exact value (dashed). The normalization chosen is a 8 = 1 and T = 0.25. The convolution 
picture always slightly overestimates the box-to-box variance (see text). 



Fig. 3. — Hierarchical mesh. The dotted lines are the buffer zones on which fine mesh initial 
seeds need to be chosen such that a seamless period-free convolution can be applied. 
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By deduction, W' R will be zero for r > R, i.e. it has compact support. The correlation kernels 
are obtained by taking the square root of the periodic correlator W PtgT ia{k) = y / ^ gr id(^) in 
Fourier space and transforming back to coordinate space. 

We will first generate initial conditions on the whole coarse grid with Wr. We linearly 
project these values onto the fine grid. This we will call 5f nc (x). We then Fourier Transform 
all the white noise conditions <5o(k) into coordinate space on the coarse grid #o( x )> an d 
also project these seeds onto the extended fine grid which we call <5Q ne (x). Let N be the 
refinement factor. The fine grid is now constant on cubicles with N 3 elements. For each 
cubicle, we will add random white noise on each grid cell of variance N 3 times the original 
white noise variance, subject to the constraint that the average must be the coarse cell 
value. Call the white noise field 5f ne . To implement the constraint, we subtract from each 
cubicle its mean and add the coarse cell value. Call the resultant density field 5| nc . We now 
convolve this field with W' R and add it to <5f ne . Since W' R is only nonzero for r < R, we 
only need to add a buffer zone of width R to the fine grid to obtain a full fine grid power 
spectrum. The actual choice of R depends on how smoothly one wishes to interpolate from 
coarse grid power to fine grid power. A typical choice of 2 grid cells should interpolate 
reasonably smoothly. 

In the language of multigrid (Press et al. 1992), we can summarize this procedure 
in terms of the prolongation operator V, the restriction operator 1Z and the convolution 
operator *. V maps from the coarse grid to the fine grid by using the nearest grid point. 
TZ does the converse map by averaging over the contained subdomain. We define the 
constrained fine grid as a white noise random field with variance < (5f nc ) 2 >= iV 3 < 5q >, 
and 5 = 1Z5f nc . Using the new operators, we can implement the constrained subpower as 
<5| ne = p<5 + 5| nc — VTZS^ nc . This gives us a compact prescription for the new fine grid 
power: 

£finc = <p WR * 5Q + W ' R * ( Wo + £finc _ pT^fine) (4) 

The same procedure can be used to embed a high resolution small volume simulation 
inside a larger grid. The reverse of the procedure above will be used. We build the large 
grid white noise field from the small grid using the restriction operator Sq = TZSq 110 in the 
region of overlap, and fill the rest of the volume with new random numbers. These are then 
convolved with the appropriate window W p . If the power on the grid scale of the small 
volume simulation was moderately small, the two grids will be consistent with each other 
in the central regions of overlap. The converse operation corresponds to cutting out a small 
volume from a larger one, and simulating the small volume at high resolution with periodic 
boundary conditions. This is straightforward if we perform the cutting on the white noise 
grid 5o( x )j and then apply the periodic convolution on the smaller white noise grid. 
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5. Conclusion 

Numerical simulations introduce unphysical truncation errors on two scales: the grid 
scale and the box length. At the grid scale, errors are dominated by the computation of 
forces and field evolution. At the box scale, currently simulations are often limited by the 
discretization of the power spectrum. For many problems in cosmology, this limit can be 
improved when described in the picture of Press-Schechter structure formation. 

We have shown how initial conditions should be generated for cosmological simulations 
utilizing small periodic box volumes. A box is considered small if the variance on the box 
scale is not small, typically greater than 10%. For CDM models, this means any box smaller 
than ~ 50/i -1 Mpc. By discretizing the correlator instead of the power spectrum, we obtain 
more accurate estimates of tophat variances and correlators, which are often the primary 
target of study in N-body simulations. 

We have also applied this approach to decompose the correlation kernel into a long 
range and short range part, allowing us to add small scale power for hierarchical subgrids 
without imposing periodicity on the subgrid subpower. This procedure generalizes to 
embedding small periodic simulations inside of larger ones and vice versa. 

I would like to thank Ben Bromley, Bill Press and Uros Seljak for useful discussions. 
This work was supported by the Harvard Society of Fellows. 
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